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Abstract 

This  survey  paper  assesses  the  status  of  compressible  Euler  and  Navier-Stokes  solvers  on  unstruc¬ 
tured  grids.  Different  spatial  and  temporal  discretization  options  for  steady  and  unsteady  flows  are 
discussed.  The  integration  of  these  components  into  an  overall  framework  to  solve  practical  problems 
is  addressed.  Issues  such  as  grid  adaptation,  higher  order  methods,  hybrid  discretizations  and  parallel 
computing  are  briefly  discussed.  Finally,  some  outstanding  issues  and  future  research  directions  are 
presented. 
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1  Introduction 


Computational  Fluid  Dynamics  (CFD)  has  evolved  rapidly  as  a  discipline  and  is  increasingly  being  used 
to  complement  the  wind  tunnel,  especially  in  preliminary  design  [50].  Based  largely  on  the  mathemat¬ 
ical  foundations  laid  among  others  by  Lax  [88]  and  Godunov  [52],  the  field  has  come  into  its  own  in 
the  last  decade.  Great  advances  have  been  made  in  the  areas  of  spatial  discretization,  grid  generation 
and  solution  strategies.  Tremendous  advances  in  computer  architecture  and  networking  speeds  have 
contributed  to  the  field  significantly  as  well.  By  the  mid  80’s,  many  different  groups  around  the  world 
w'ere  able  to  compute  three-dimensional  flows  over  simple  realistic  aerodynamic  configurations.  The 
grids  employed  were  of  the  body-fitted  or  structured  grid  type.  One-dimensional  models  were  extended 
to  deal  with  multiple  dimensions  in  a  natural  way  because  of  the  structure  by  using  the  so-called  gener¬ 
alized  coordinates  [161].  However,  the  task  of  generating  structured  grids  about  complex  configurations 
presented  a  serious  challenge.  The  widely-used  multi-block  structured  grid  approach  solves  this  problem 
by  tessellating  the  domain  between  the  body  and  the  far- field  into  simple  logically  rectangular  blocks, 
so  that  structured  grids  can  be  generated  easily  within  each  block  [17.3,  81].  The  automation  of  the 
blocking  and  the  grid  generation  process  are  difficult  tasks  that  are  continually  being  refined.  Another 
powerful  approach  uses  overlapping  or  chimera  grids  [19,  110].  Here,  structured  grids  generated  about 
the  different  components,  are  allowed  to  overlap.  Automation  of  blocking,  grid  generation  and  the 
preprocessing  required  for  deriving  the  interpolation  operators  are  continually  being  improved. 

The  desire  to  compute  flows  over  complex  configurations  also  spawned  a  surge  of  activity  in  the 
area  of  unstructured  grids.  The  term  “unstructured  grids”  will  be  used  primarily  in  this  paper  to  mean 
grids  composed  of  simplices.  which  are  triangles  in  two  dimensions  and  tetrahedra  in  three  dimensions. 
Unstructured  grids  have  always  been  used  in  finite  element  circles,  but  have  become  popular  in  the 
finite  volume  community  only  fairly  recently.  They  provide  flexibility  for  tessellating  about  complex 
geometries  and  for  adapting  to  flow  features,  such  as  shocks  and  boundary  layers.  An  underlying 
premise  is  that  unstructured  grid  generation  is  far  more  automatable  than  are  the  tasks  associated  with 
multi-block  structured  grid  generation. 

It  should  be  mentioned  that  the  flow  solver  only  constitutes  a  part  of  the  overall  solution  method¬ 
ology.  Other  tasks,  such  as  grid  generation,  interfacing  with  geometry  packages,  visualization  and 
post-processing  are  equally  important,  and  demand  considerable  skills  and  resources.  The  outline  of 
this  survey  paper  on  unstructured  grid  flow  solvers  is  as  follows.  Section  2  presents  the  governing 
equations  in  integral  form.  Section  3  reviews  some  popular  finite  volume  spatial  discretization  schemes. 
Section  4  reviews  the  finite  element  approach  for  spatial  discretization.  Section  5  reviews  the  turbulence 
models  in  vogue  for  unstructured  grid  computations.  Section  6  presents  the  various  alternatives  avail¬ 
able  for  time-discretization  for  steady  and  unsteady  problems.  Section  7  addresses  the  important  topic 
of  grid  adaptation.  Sections  8  and  9  examine  higher  order  accurate  schemes  and  hybrid  discretizations. 
Section  10  briefly  examines  parallel  computing  issues.  The  paper  concludes  with  a  listing  of  issues  to  be 
resolved  and  some  future  research  directions.  In  addition  to  the  papers  cited  in  this  article,  the  reader  is 
encouraged  to  read,  especially  regarding  spatial  discretization,  the  AGARD  Report  [1]  on  unstructured 
grid  methods  for  advection  dominated  flows. 

2  Governing  equations 

The  equations  governing  compressible  fluid  flow  in  integral  form  for  a  control  volume  V(t)  with  boundary 
S{t)  are  given  by 

^  [  Wdv  +  i  [F{W,  n,  s)  -  G{W,  VVF,  n)]da  =  0,  (1) 
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In  the  formulas  given  above  p  is  the  density,  V  is  the  velocity  vector  with  Cartesian  components  e  is 
the  specific  total  energy,  n  is  the  outward  unit  normal  vector  of  the  boundary  S{  t)  and  s  is  the  velocity 
vector  of  the  boundary.  Also,  p  is  the  molecular  viscosity,  A  is  the  bulk  viscosity  related  to  p  by  Stokes’ 
hypothesis,  A  =  -2/.3/i,  7  is  the  identity  tensor,  T  is  the  stress  tensor  and  D  is  the  deformation  tensor 
given  by 

D.,  =  +  Vi,i)  (2) 


where  Vi  j  denotes  the  partial  derivative  of  the  z’th  component  of  V  with  respect  to  the  Cartesian 
coordinate  x^,  i.e.  =  |^-  0  stands  for  the  divergence  of  V  given  by  V^i  with  the  usual  summation 
convention,  q  is  the  heat  flux  given  by  Fourier’s  law 


q  =  -KVT\ 


(3) 


where  K  is  the  thermal  conductivity  of  the  fluid  and  T  is  the  temperature.  These  equations  are 
augmented  by  the  equation  of  state,  which  for  a  perfect  gas  is  given  by 

p  =  h-l)(pe-tp|Vp)  (4) 

Eqn.  (1)  represents  the  conservation  laws  for  the  mass,  momentum  (the  Navier-Stokes  equations)  and 
energy.  It  holds  for  any  volume  and  in  particular,  holds  for  a  specific  volume  associated  with  each  grid 
point,  termed  the  control  volume. 


3  Finite  volume  spatial  discretization 

It  is  assumed  that  a  grid  about  the  geometry  of  interest  has  been  generated  by  some  suitable  method. 
The  various  grid  generation  techniques  in  use  are  reviewed  in  the  survey  papers  by  Thompson  and 
Weatherill  [166]  and  Mavriplis  [106].  On  a  given  grid,  one  has  at  least  two  choices  as  to  where  to  locate 
the  variables,  giving  rise  to  the  cell- vertex  and  the  cell-centered  approaches.  In  the  cell-vertex  approach, 
the  variables  are  stored  at  the  vertices  of  the  grid,  whereas  in  the  cell-centered  approach  they  are  stored 
at  the  centroids  of  the  cells.  There  is  yet  another  approach  that  stores  only  the  averages  associated 
with  control  volumes.  If  the  scheme  possesses  second  order  spatial  accuracy  or  less,  this  approach  is  no 
different  than  a  cell-centered  approach;  the  higher  order  scheme  is  covered  briefly  in  Section  8. 

The  concept  of  using  arbitrary  control  volumes  to  solve  numerically  the  conservation  laws  was  es¬ 
tablished  by  the  late  70’s,  at  least  in  theory  [89].  Jameson  and  Mavriplis  [73]  reported  some  of  the 
earliest  results  from  solving  the  two-dimensional  Euler  equations  on  regular  triangular  grids  that  were 
obtained  by  subdividing  quadrilateral  grids.  In  a  cell-centered  setting,  they  extended  much  of  what  was 
established  for  finite  volume  schemes  on  structured  grids  [74,  67]  to  triangular  grids.  This  included  the 
constructions  of  central-difference-like  approximation  for  the  convective  terms  and  a  blend  of  dissipative 
terms  to  suppress  odd-even  decoupling  and  to  capture  discontinuities,  and  the  incorporation  of  multi¬ 
grid  ideas.  Solutions  and  convergence  rates  of  comparable  quality  to  structured  grids  were  obtained. 
Second  order  accuracy  was  demonstrated  by  using  a  sequence  of  uniform  grids.  In  1986,  Jameson  et 
al.  [72]  presented  their  paper  dealing  with  inviscid  transonic  flow  over  a  complete  aircraft.  The  pa¬ 
per’s  contributions  included  grid  generation  for  complex  geometries  using  the  Delaunay  triangulation 
approach  and  the  development  of  a  cell- vertex  flow  solver.  The  control  volume  for  each  vertex  in  the 
tetrahedral  grid  was  taken  to  be  the  union  of  all  tetrahedra  sharing  that  vertex.  They  also  showed 
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the  correspondence  between  this  and  a  Galerkin  finite  element  procedure  with  piecewise-linear  basis 
functions.  The  nondissipative  approximation  was  augmented  with  dissipative  terms.  In  addition,  the 
paper  also  presented  a  first  discussion  on  a  scheme  with  positive  coefficients  for  triangles,  although  the 
proposed  scheme  was  only  first  order  accurate.  Figure  1  shows  the  surface  Mach  contours  for  inviscid 
transonic  flow  over  a  Boeing  747-200  taken  from  the  paper  of  Jameson  et  al.  [72].  The  flow  conditions 
are  Moo  =  0.8  and  a  -  2.73°.  The  mesh  used  for  the  calculation  contained  12,038  nodes  (-57914  tetra- 
hedra)  and  the  calculation  took  about  25  minutes  on  a  Cray  XMP  computer.  Since  this  seminal  effort, 
tremendous  advances  have  been  made  in  grid  generation  and  vectorization  of  the  flow  solver  and  are 
outlined  in  the  sequel  [71]. 

In  1973,  Boris  and  Book  [25]  introduced  the  Flux  Corrected  Transport  (FCT)  scheme  that  converted 
a  first  order  accurate  monotone  scheme  to  a  higher  order  scheme  by  adding  limited  amounts  of  anti- 
diffusive  flux  to  prevent  spurious  oscillations.  Harten  [54]  developed  a  mathematical  formulation  that 
gave  rise  to  the  Total  Variation  Diminishing  (TVD)  schemes  which  also  added  limited  amounts  of  anti- 
diffusive  fluxes  to  a  first  order  scheme  to  achieve  monotonicity.  Van  Leer  [169]  independently  devised  the 
concept  of  limiting  when  designing  monotonicity  preserving  schemes  for  conservation  laws.  He  devised 
the  Monotonic  Upstream-centered  Scheme  for  Conservation  Laws  (MUSCL)  schemes  which  relied  on  a 
piecewise-polynomial  reconstruction  procedure  that  enforced  monotonicity  principles  by  using  nonlinear 
functions  called  limiters.  The  jumps  at  the  interfaces  were  resolved  by  using  a  Riemann  solver,  much  like 
Godunov’s  scheme.  Since  the  Riemann  problem  for  the  Euler  equations  requires  an  iterative  procedure 
and  is  expensive,  the  quest  was  on  for  approximate  Riemann  solvers.  Roe  [142,  143]  and  Osher  [120] 
constructed  two  such  schemes  that  proved  to  be  particularly  good  in  concert  with  the  MUSCL  approach 
on  structured  grids.  For  an  excellent  discussion  on  the  development  of  upwind  schemes,  see  the  text  by 
Leveque  [91].  By  the  mid  80’s,  MUSCL  schemes  were  being  used  by  a  number  of  groups  to  compute 
aerodynamic  flows  on  structured  grids  [164,  160].  These  schemes  were  gradually  being  extended  to  deal 
with  unstructured  grids. 

Desideri  and  Dervieux  [40]  devised  cell-vertex  finite  volume  schemes  for  unstructured  grids  using 
MUSCL  ideas.  Given  pointwise  values  at  the  ceU  vertices  of  a  triangulation,  they  employed  a  recon¬ 
struction  procedure  that  made  use  of  gradients  in  neighboring  triangles.  The  flux  was  constructed  using 
Osher’s  approximate  Riemann  solver.  Limiters  were  used  in  a  one- dimensional  fashion  but  a  multi¬ 
dimensional  monotonicity  principle  was  not  satisfied.  Nevertheless,  respectable  results  were  obtained 
for  many  inviscid  computations.  Lohner  et  al.  [95]  tested  a  FEM-FCT  scheme  for  Euler  and  Navier- 
Stokes  equations.  Fezoui  and  Stoufflet  [4-5]  proposed  and  tested  a  class  of  implicit  upwind  schemes 
that  utilized  various  upwind  approximations.  Whitaker  et  al.  [185]  constructed  a  similar  scheme  using 
Roe’s  flux-difference  splitting  scheme  and  obtained  results  for  many  transonic  and  supersonic  flows. 
Venkatakrishnan  and  Barth  [176]  constructed  an  upwind  scheme  for  cell-centered  triangular  grids  that 
also  employed  MUSCL  ideas  as  did  Batina  [16],  Knight  [84]  and  Frink  [47].  However,  none  of  these 
efforts  reaUy  guaranteed  the  absence  of  oscillations  in  the  multi-dimensional  case. 

Barth  and  Jespersen  [15]  made  a  radical  departure  from  one- dimensional  thinking  for  satisfying 
monotonicity  principles.  They  enunciated  a  monotonicity  principle  in  multiple  dimensions  similar  to 
that  employed  by  Van  Leer  [169]  and  Spekreijse  for  structured  grids  [160],  namely  that  the  reconstructed 
distribution  in  the  control  volume  be  bounded  everywhere  by  the  values  of  the  neighbors  (including 
the  vertex  representing  the  control  volume)  and  satisfied  the  principle  by  constructing  a  truly  multi¬ 
dimensional  limiter.  They  employed  Roe’s  approximate  Riemann  solver  for  the  evolution  phase.  As 
an  ultimate  test,  they  computed  oscillation-free  solutions  for  transonic  flow  over  an  airfoil  on  a  highly 
irregular  mesh.  They  obtained  solutions  of  comparable  accuracy  with  both  cell-centered  and  cell- vertex 
schemes.  This  MUSCL  approach  has  been  adopted  since  by  other  researchers  [5,  7,  51].  The  multi¬ 
dimensional  limiter  in  [15]  may  be  thought  of  as  a  generalization  of  the  min-mod  limiter  and  as  such, 
leads  to  convergence  difficulties.  Venkatakrishnan  [174]  has  analyzed  this  problem  and  has  proposed 
modifications  that  ameliorate  the  situation  at  the  expense  of  monotonicity.  Aftosmis  et  al.  [-3]  have 
found  that  the  modifications  proposed  in  [174]  significantly  improve  the  convergence  as  well  as  solution 
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accuracy  in  test  problems.  Shu  [15-3]  has  proposed  the  Total  Variation  Bounded  (TVB)  schemes  which 
weaken  the  TV^D  property  by  allowing  for  small-scale  oscillations  although  the  motivation  here  is  to 
avoid  clipping  smooth  extrema.  The  effect  of  TVB  modification  on  convergence  has  not  been  studied 
yet . 

Frink  et  al.  [47,  49]  developed  an  upwind  cell-centered  three-dimensional  flow  solver.  They  employed 
a  weighted  averaging  procedure  that  interpola.ted  variables  from  the  centers  of  tetrahedra  to  the  vertices 
and  used  these  vertex  values  to  compute  the  gradients  within  each  tetrahedron.  These  gradients  were 
used  to  interpolate  the  variables  to  the  centers  of  the  faces  of  the  tetrahedra.  This  was  followed  by  the 
use  of  Roe’s  approximate  Riemann  solver.  The  reconstruction  procedure  was  linear  and  therefore,  not 
monotonicity  preserving.  Nevertheless,  the  averaging  procedure  seemed  to  introduce  enough  dissipation 
that  good  transonic  flow  solutions  were  obtained  without  using  limiters.  The  reconstruction  procedure 
has  since  been  modified  to  be  linearity-preserving  [48]  based  on  the  work  of  Holmes  and  ConneU  [61]. 
Laminar  viscous  capability  has  been  been  added  to  this  code.  Figure  2a  depicts  the  surface  grid  and 
“oil  flow”  pattern  for  laminar  flow  over  a  delta  wing  at  Moo  =  0.3,  a  =  20.5°,  Rei  =  0.9  x  10®. 
The  flow  has  been  computed  by  Frink  using  the  methodology  described  in  [48].  The  grid  containing 
730,454  tetrahedral  cells  was  generated  by  the  Advancing-Layers  Method  [132]  which  produces  thin- 
layer  tetrahedral  grids  suitable  for  computing  viscous  flows.  Evidence  of  primary,  secondary  and  tertiary 
vortices  may  be  seen  in  Figure  2a..  Figure  2b  shows  the  pressure  profiles  at  four  chord  stations  and  the 
comparisons  with  the  results  from  the  structured  grid  code  CFL3D  [164]  and  the  experimental  data  of 
Hummel  [66]. 

Barth  [11,  12]  presented  extensions  of  the  schemes  of  [15]  to  three-dimensional  inviscid  and  viscous 
flow  computations.  The  main  contributions  of  these  papers  included  the  discretization  of  the  viscous 
terms  using  a  finite  element  procedure  and  the  use  of  edge-based  data  structures  for  discretizing  the 
inviscid  and  viscous  terms.  In  addition,  least-squares  and  data- weighted  procedures  for  the  construction 
of  gradients  were  outlined.  Mavriplis  [102,  105,  101]  developed  an  explicit  vertex-based  finite  element 
multigrid  scheme  for  the  two-  and  three-dimensional  Euler  and  Navier-Stokes  equations.  He  also  pre¬ 
sented  edge-based  data  structures.  In  [11, 105,  98,  97]  a  case  is  made  for  using  cell- vertex  approximations 
and  edge-based  data  structures  in  three  dimensions  by  examining  the  computational  complexity  and 
storage. 

The  issue  of  cell-vertex  vs  cell-centered  approximations  is  still  an  open  one,  particularly  in  three 
dimensions.  In  two  dimensions  the  ratio  of  the  number  of  cells  to  the  number  of  vertices  is  2  whereas 
in  three  dimensions  this  ratio  could  be  arbitrarily  large,  although  it  is  typically  between  5  and  6  for 
nice  tetrahedralizations.  In  the  case  of  tetrahedral  grids,  the  flux  computations  in  a  cell-vertex  scheme 
can  be  cast  as  loops  over  edges  whereas  in  a  cell-centered  scheme  they  are  loops  over  triangular  faces. 
The  ratio  of  the  number  of  faces  to  the  number  of  edges  in  a  typical  tetrahedral  grid  is  roughly  2. 
Therefore,  it  is  true,  as  argued  in  [11,  105]  that  on  a  given  grid,  the  cell-centered  approximation  incurs 
considerably  more  computational  effort  and  memory  requirements.  However,  there  is  some  evidence 
that  on  a  given  grid,  the  solution  quality  of  a  cell-centered  scheme  is  superior  to  that  of  a  cell-vertex 
scheme  [134].  This  is  more  than  likely  due  to  the  fact  that  for  triangular/tetrahedral  grids  the  control 
volumes  in  a  cell-centered  scheme  are  smaller  than  those  in  a  cell- vertex  scheme.  It  is  not  clear  whether 
the  cell- vertex  scheme  requires  a  grid  that  has  as  many  vertices  as  the  number  of  tetrahedra  employed 
by  a  cell-centered  scheme  to  achieve  the  same  level  of  accuracy.  If  this  is  the  case,  the  cell-centered 
scheme  wiU  prevail  because  it  requires  much  smaller  grids  to  be  generated.  However,  it  appears  that  a 
cell-vertex  scheme  is  better  suited  to  computing  the  viscous  fluxes,  especially  when  the  triangulations 
become  highly  irregular.  As  discussed  later  in  this  section,  the  viscous  terms  in  a  cell-vertex  scheme 
are  typically  discretized  using  a  Galerkin  finite  element  approach,  whereas  a  finite- volume  viewpoint 
is  adopted  in  a  cell-centered  scheme.  It  has  been  shown  [11,  90]  that  a  finite  element  discretization 
of  a  diffusion  operator  with  linear  basis  functions  obeys  the  discrete  maximum  principle  if  a  Delaunay 
triangulation  is  used.  This  result  is  true  only  in  two  dimensions  and  does  not  hold  in  three  dimensions. 
On  the  other  hand,  the  finite  volume  discretization  of  the  diffusion  operator  in  a  cell-centered  setting 
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does  not  guarantee  a  discrete  maximum  principle  in  two  or  three  dimensions.  It  is  interesting  to  note 
that  the  centroidal  control  volumes  that  are  used  with  cell-vertex  schemes  are  not  necessarily  convex 
and  can  assume  rather  odd  shapes,  especially  in  the  case  of  highly  stretched  grids.  One  alternative  is 
to  use  the  Voronoi  cells,  which  are  guaranteed  to  be  convex  [86],  but  this  approach  requires  care  at  the 
boundaries. 

Significant  developments  have  also  taken  place  in  the  area  of  positive  schemes.  These  require  that 
the  semi-discrete  scheme  be  expressible  in  the  form 

where  the  Cij’s  are  non-negative,  and  Ni  denotes  the  set  of  neighbors  of  i.  As  observed  in  [72],  if  the 
C'ij's  are  constants,  the  scheme  is  linear  and  hence  only  first  order  accurate.  Higher  order  versions 
of  such  schemes  have  been  developed  [42,  .33,  69,  70].  These  may  be  regarded  as  generalizations  of 
TVD  schemes  to  conservation  laws  in  multi-dimensions.  Durlofsky  et  al.  [42]  have  developed  a  cell- 
centered  scheme  that  makes  use  of  gradients  in  auxiliary  triangles  similar  to  the  scheme  of  [176],  but 
incorporates  a  multi- dimensional  limiter  to  prevent  oscillations.  Both  the  schemes  have  the  drawback 
of  requiring  nice  triangulations,  since  the  auxiliary  gradients  could  be  ill-defined  on  highly  nonuniform 
meshes.  Jameson  [69]  has  developed  a  theory  for  Local  Extremum  Diminishing  (LED)  and  Essentially 
Local  Extremum  Diminishing  (ELED)  schemes  leading  to  Symmetric  and  Upstream  Limited  Positive 
(SLIP  and  USLIP)  schemes.  He  has  also  presented  an  analysis  for  schemes  to  resolve  shocks  with  one 
interior  point  [70].  A  simple  scheme  that  meets  the  conditions  is  termed  a  Convective  Upwind  and  Split 
Pressure  (CUSP)  scheme,  which  is  closely  related  to  the  AUSM  scheme  of  Steffen  and  Liou  [92].  The 
schemes  have  been  tested  on  structured  grids  for  inviscid  and  viscous  flows  [162]. 

A  different  avenue  to  realizing  upwinding  on  unstructured  grids  is  provided  by  the  fluctuation¬ 
splitting  schemes  [39,  121].  In  multi-dimensions,  the  weak  link  in  the  MUSCL  formulation  is  the  use  of 
a  directionally-split  Riemann  solver.  This  approach  is  routinely  adopted,  but  can  misinterpret  features 
not  aligned  with  the  grid.  To  overcome  this  problem,  various  multi-dimensional  Riemann  solvers  have 
been  proposed  in  the  literature  [38,  12.5,  148,  146],  but  have  not  gained  widespread  acceptance.  Rather 
than  adopt  the  approach  of  reconstruction  followed  by  an  approximate  Riemann  solver  to  bring  in 
the  upwinding,  fluctuation-splitting  schemes  consider  the  average  time  evolution  of  a  complete  cell  (a 
triangle  or  a  tetrahedron)  with  the  unknowns  located  at  its  vertices,  and  then  update  the  values  at 
the  vertices  by  the  effect  of  linear  wave  solutions  evolving  the  piecewise  linear  data  over  a  cell.  If  the 
distribution  formulas  are  unbiased,  this  leads  to  the  well-known  Ni’s  scheme  [119],  a  centered  scheme 
which  requires  the  addition  of  dissipative  terms.  Several  uwpind  fluctuation-splitting  schemes  for  scalar 
advection  equation  have  been  developed;  these  are  compared  in  [121].  Extension  of  fluctuation-splitting 
schemes  to  systems  of  equations  relies  on  the  decomposition  of  flux  divergence  into  scalar  waves.  The 
decomposition  is  not  unique,  and  several  strategies  such  as  characteristic  decomposition  [38]  and  simple 
wave  models  [144,  125]  have  been  investigated.  Some  of  the  advantages  of  the  fluctuation-splitting 
approach  are  that  they  result  in  compact  stencils  and  that  they  do  not  use  ad-hoc  dimensionaUy-split 
Riemann  solvers.  Some  of  the  disadvantages  include  the  fact  that  the  scheme  is  only  second  order 
accurate  at  steady  state,  and  that  they  bind  the  user  to  using  a  decomposition  technique.  Recently,  a 
positive  linearity-preserving  fluctuation-splitting  scheme  has  been  constructed  for  the  multi-dimensional 
scalar  advection  equation  using  limiters  [155].  It  differs  from  the  other  approaches  in  that  a  purely 
algebraic  (not  geometric)  viewpoint  is  adopted.  This  approach  has  enabled  Sidilkover  [154]  to  devise  a 
fluctuation-splitting  scheme  for  the  Euler  equations  that  may  be  viewed  as  a  genuine  multi-dimensional 
extension  of  Roe’s  scheme  without  having  to  resort  to  decomposition  into  scalar  waves.  A  unique 
property  of  this  scheme  is  that  a  nonlinear  Gauss  Seidel  procedure  can  be  used  as  a  relaxation  scheme. 

In  finite  volume  schemes,  a  spatial  (possibly  monotonicity-preserving)  approximation  for  the  inviscid 
terms  is  combined  with  a  centered  approximation  for  the  viscous  fluxes.  A  basis  for  this  “operator¬ 
splitting”  approach  is  furnished  by  Roe[14.5]  who  showed  by  means  of  elegant  analysis  in  one  dimension 
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how  a  monotonicity-preserving  discretization  for  the  inviscid  terms  when  combined  with  a  centered 
approximation  for  the  viscous  terms  suppresses  spurious  modes  and  also  yields  a  nonoscillatory  solution 
for  all  cell  Reynolds  numbers.  On  unstructured  grids,  the  viscous  fluxes  are  discretized  by  using  either 
finite  volume  or  finite  element  approaches.  In  the  case  of  cell-centered  discretizations  Frink  et  al.  [48] 
compute  the  first  derivatives  at  the  vertices  of  the  triangulation,  which  are  then  averaged  to  obtain  the 
viscous  fluxes  at  the  faces.  In  the  case  of  cell-vertex  schemes,  the  viscous  terms  are  usually  derived  by 
adopting  a  Galerkin  finite  element  approach  [147,  11,  101],  although  there  is  an  equivalent  finite  volume 
interpretation.  This  formulation  results  in  a  compact  stencil  involving  only  the  nearest  neighbors,  and 
also  exactly  reproduces  the  gradient  of  a  linear  function.  A  simpler  formulation  has  been  proposed  in 
[98].  Here,  the  viscous  fluxes  are  computed  by  averaging  the  cell-vertex  gradients,  which  are  available 
from  the  reconstruction  procedure  for  the  inviscid  terms.  However,  this  method  results  in  a  larger 
stencil  for  the  viscous  terms  compared  to  the  finite  element  discretization  and  could  be  more  diffusive. 

4  Finite  element  discretizations 

Finite  element  method  (FEM)  merits  a  separate  discussion  since  its  evolution  in  CFD  has  taken  place  in 
parallel  with  finite  volume  schemes.  The  rigorous  treatment  adopted  in  FEM  has  important  implications 
for  global  error  estimation  and  grid  adaptation.  Many  researchers  have  shown  the  equivalence  between 
finite  volume  and  finite  element  approaches  when  piecewise-linear  approximations  are  employed  for  the 
solution  vector  and  the  fiuxes  [147,  72,  11,  150]. 

The  most  powerful  method  in  finite  elements  is  the  Galerkin  method.  Here,  the  solution  is  first 
expanded  in  a  set  of  basis  functions  and  the  residual  is  made  orthogonal  to  a  set  of  test  functions.  When 
the  basis  and  the  test  functions  are  the  same,  the  method  is  termed  a  Galerkin  method;  otherwise  it 
is  called  a  Petrov- Galerkin  method.  The  standard  Galerkin  method  leads  to  a  centered  scheme  and  is 
unconditionally  unstable  for  hyperbolic  problems  when  combined  with  forward  Euler  discretization  in 
time.  Therefore,  artificial  viscosity  has  to  be  added  in  some  form  to  stabilize  the  procedure. 

In  contrast  with  the  finite  volume  methods,  finite  element  practitioners  have  always  preferred  cell- 
vertex  schemes  since  the  global  function  is  usually  expressed  as  a  summation  of  trial  functions  multiplied 
by  the  values  at  the  vertices.  These  trial  functions  are  typically  assumed  to  be  1  at  a  vertex  and  zero  at 
ail  other  vertices.  Bristeau  et  al.  [28]  computed  many  low  Reynolds  number  compressible  flows  using  a 
finite  element  method.  Angrand  and  Dervieux  [8]  investigated  several  first  and  second  order  accurate 
explicit  schemes.  They  added  Lapidus-type  artificial  viscosity  to  a  Galerkin  method.  Donea  [41]  derived 
the  important  Taylor- Galerkin  family  of  schemes  for  the  linear  advection  equation.  Adopting  an  FEM 
approach,  he  showed  how  a  Galerkin  scheme  (a  centered  scheme)  could  be  stabilized  by  using  a  Taylor- 
series  expansion  for  the  time  derivative  ||,  similar  to  the  procedure  used  to  derive  the  Lax-Wendroff 
scheme.  He  demonstrated  that  the  resulting  schemes  had  good  phase  error  and  dissipation  properties 
and  that  they  could  be  easily  extended  to  multi-dimensions.  Lohner  et  al.  [96]  developed  and  tested  a 
two-step  Taylor- Galerkin  finite  element  method  for  the  solution  of  the  Euler  equations.  Since  the  use 
of  artificial  viscosity  spread  the  discontinuities  over  several  cells,  they  also  employed  an  adaptive  grid 
algorithm  using  local  refinement. 

While  there  are  finite  difference/ volume  equivalents  to  most  of  the  finite  element  techniques  des¬ 
cribed  in  the  previous  paragraph,  a  very  important  class  of  schemes  has  been  derived  in  the  finite  element 
community  that  has  no  obvious  counterpart  in  the  finite  difference/ volume  approach.  Included  in  this 
class  are  the  Streamwise  Upwind  Petrov  Galerkin  (SUPG)  or  the  Streamwise  Diffusion  (SD),  and  the 
Galerkin  Least  Squares  methods.  SUPG,  devised  originally  by  Hughes  and  Brooks  [63]  for  the  steady 
scalar  advection- diffusion  equation,  and  Galerkin  least  squares  [64]  methods  have  subsequently  been  ex¬ 
tended  to  deal  with  the  compressible  and  incompressible  Navier-Stokes  equations  [62].  To  the  Galerkin 
discretization,  they  add  terms  proportional  to  the  residual  (unsteady  residual  for  time-dependent  prob¬ 
lems)  to  introduce  dissipation.  Unlike  conventional  artificial  viscosity,  this  form  of  artificial  viscosity 
does  not  compromise  the  order  of  accuracy  of  the  scheme  since  the  exact  solution  yields  zero  dissipation. 
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The  global  error  in  L2  norm  can  been  shown  to  vary  as  for  the  advection-dominated  case  and  as 

for  the  diffusion- dominated  case,  where  p  is  the  degree  of  piecewise-polynomial  basis  function  [62]. 
Another  method  that  has  been  extensively  analyzed  is  the  discontinuous  Galerkin  method  [77,  .33,  22], 
•Johnson  et  al.  [78]  have  analyzed  the  SUPG  scheme  and  have  extended  it  to  deal  with  time-dependent 
problems  by  adopting  a  discontinuous  Galerkin  procedure  in  time.  A  review  and  analysis  of  these  and 
other  finite  element  methods  may  be  found  in  the  text  by  .Johnson  [76].  All  these  methods  are  hnear  and 
are  not  suitable  for  capturing  discontinuities.  Hughes  and  Mallet  [6.5]  have  proposed  a  shock-capturing 
operator  that  locally  reduces  the  order  of  accuracy  near  discontinuities.  Three-dimensional  viscous  flow 
over  the  canopy  of  a  Hermes  computed  by  Chalot  et  al.  [32]  is  presented  in  Figures  3(a-b).  The  flow 
conditions  corresponding  to  reentry  are  altitude  of  60  km.  Moo  —  20,  angle  of  attack  a  =  30°  and  Re/m 
of  120,890.  The  laminar  solution  was  computed  using  equilibrium  real  gas  hypothesis  by  the  SUPG 
finite  element  method  by  an  iterative  implicit  method  that  employs  a  block-diagonal  preconditioned 
GMRES  procedure.  The  mesh  contains  nearly  1  million  tetrahedra.  Figure  .3a  displays  the  surface  grid 
and  Figure  3b  depicts  the  skin  friction  lines,  illustrating  the  complex  flow  structure. 


5  Turbulence  modeling 

In  order  to  compute  practical  flows  at  large  Reynolds  numbers,  turbulence  has  to  be  modeled.  Some  of 
the  most  popular  turbulence  models  in  use  in  structured  grids  such  as  the  Baldwin-Lomax  [10]  become 
quite  difficult  to  implement  on  unstructured  grids  because  of  their  nonlocal  nature,  although  this  has 
been  done  [147,  104].  The  trend  is  away  from  algebraic  models  and  towards  simple  one  and  two  equation 
field  models,  k  -  e  turbulence  models  have  been  used  with  unstructured  grids  either  by  integrating  the 
model  to  the  wall  [168]  or  by  using  a  two-layer  model  [116,  32],  where  the  two  equation  k  -  e  model  is 
replaced  by  a  one  equation  model  near  the  wall  region.  Marcum  and  Agarwal  [100]  have  implemented 
two  versions  of  k  -  e  models  in  an  unstructured  grid  code  to  compute  turbulent  flows  over  an  ONERA 
M6  wing.  They  tested  a  low  Reynolds  number  model  that  is  integrated  to  the  wall  and  a  high  Reynolds 
number  model  that  is  combined  with  wall-functions.  Two  fairly  new  one-equation  models  have  become 
very  popular,  particularly  for  unstructured  grid  applications.  These  are  the  Baldwin-Barth  [9]  and  the 
Spalart-AUmaras  [159]  models.  The  Baldwin-Barth  model  is  derived  from  a  simplified  form  of  the  k-e 
turbulence  model.  The  Spalart-AUmaras  model,  on  the  other  hand,  is  derived  “from  scratch”  through 
dimensional  analysis  and  intuition.  Both  the  models  have  been  tested  extensively  and  are  being  used 
routinely.  The  cell-vertex  finite  volume  scheme  with  a  Galerkin  procedure  for  computing  the  viscous 
terms  has  been  used  to  compute  viscous  flows  about  high-lift  configurations  in  conjunction  with  these 
turbulence  models  [6,  168].  In  Figure  4,  the  surface  pressure  distributions  and  velocity  profiles  at  three 
stations  on  the  flap  are  displayed.  The  results  are  taken  from  the  paper  of  Anderson  and  Bonhaus  [6]. 
They  use  an  implicit  cell- vertex  MUSCL  scheme  with  the  Spalart-AUmaras  turbulence  model,  where  the 
linear  system  at  each  time  step  is  solved  by  a  Gauss  Seidel  relaxation  method.  The  results  are  compared 
with  experimental  data  at  two  Reynolds  numbers  5  x  10®  and  9  x  10®  with  =  0.2  and  angle  of 
attack,  a  =  16°.  The  surface  pressure  profiles  show  excellent  agreement  with  the  experimental  data, 
while  the  velocity  profiles  show  good  agreement  except  for  the  wake  regions  far  downstream,  where  the 
grid  resolution  is  inadequate.  These  and  other  quantitative  comparisons  with  experimental  data  [168] 
indicate  the  level  of  maturity  of  unstructured  grid  technology. 


6  Time  Discretization  and  Solution  Strategies 


6.1  Steady-state  solution  techniques 

After  discretizing  Eqn.  (1)  in  space,  the  following 
obtained: 

d{VMW) 

dt 


system  of  coupled  ordinary  differential  equations  is 
d-  RiW)  =  0.  (6) 
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Here  W  is  the  vector  of  unknowns  over  the  mesh  points  for  a  vertex-based  formulation  and  over  the 
cells  for  a  cell-based  formulation.  V  is  the  volume  of  the  polyhedral  control  volume  associated  with 
the  vertex/ceU.  In  the  case  of  a  cell-vertex  scheme  storing  pointwise  quantities  at  the  vertices,  M  is 
the  mass  matrix  which  represents  the  relationship  between  the  average  value  in  a  control  volume  and 
the  values  at  the  vertices  (the  vertex  representing  the  control  volume  and  its  nearest  neighbors).  It  is 
only  a  function  of  the  mesh  and  hence,  a  constant  matrix  for  a  static  mesh.  If  a  steady  state  solution 
is  sought,  time- accuracy  is  not  an  issue  and  M  can  be  replaced  by  the  identity  matrix.  This  technique 
known  as  mass-lumping  yields  the  following  system  of  ordinary  differential  equations  for  the  vector  of 
unknowns  W  : 

+  Ji(W)  =  0.  (7) 

Cell-centered  schemes  (up  to  second  order  accuracy)  and  schemes  dealing  strictly  with  cell-averages  (to 
any  order  of  accuracy)  do  not  yield  a  mass  matrix  and  thus  lead  to  Eqn.  (7)  for  steady  and  unsteady 
problems. 

Eqn.  7  may  be  solved  explicitly  with  linear  multi-step  methods  of  the  forms  described  in  [6y, 
170].  The  coefficients  for  these  Runge-Kutta  schemes  are  derived  by  considering  a  model  problem  and 
optimizing  the  coefficients  to  yield  a  large  CFL  number  and  good  damping  properties.  Local  time 
stepping,  enthalpy  damping  (for  Euler  equations)  and  residual  averaging  are  employed  to  accelerate 
convergence  [74].  Even  with  this  methodology,  the  convergence  to  steady  state  is  usually  unacceptably 
slow.  Therefore,  either  multigrid  methods  or  implicit  schemes  are  required  to  accelerate  the  convergence. 

6.1.1  Multigrid  methods 

The  multigrid  method  has  been  demonstrated  as  an  efficient  means  for  obtaining  steady-state  solutions 
to  the  compressible  Euler  and  Navier- Stokes  equations  on  unstructured  meshes  in  two  and  three  dimen¬ 
sions.  In  this  approach,  convergence  acceleration  is  achieved  by  time-stepping  on  successively  coarser 
meshes.  The  principle  behind  this  algorithm  is  that  the  errors  associated  with  the  high  frequencies 
are  damped  by  a  carefully  chosen  smoother  (e.g.  a  multi-stage  Runge-Kutta  scheme)  while  the  errors 
associated  with  the  low  frequencies  are  damped  on  the  coarser  grids  where  these  frequencies  manifest 
themselves  as  high  frequencies.  In  the  case  of  structured  grids,  coarse  grids  are  easily  derived  from  a 
given  fine  grid  by  omitting  alternate  grid  fines  in  each  coordinate  direction.  In  the  case  of  unstructured 
grids,  three  dilferent  approaches  can  be  adopted. 

The  first  approach  begins  with  a  coarse  mesh  definition  and  generates  finer  grids  by  refinement 
[1.31,  35].  The  advantage  is  that  the  inter-grid  operators  become  simple  because  of  the  nesting  of  grids. 
Another  advantage  is  that  this  set-up  can  be  utilized  to  advantage  in  an  adaptive  procedure,  where  the 
fine  meshes  are  formed  by  adaptively  refining  the  coarse  meshes  [35,  127].  The  principal  disadvantage  of 
this  approach  is  the  dependence  of  the  fine  grid  distribution  on  the  coarse  levels.  The  second  approach 
uses  non-nested  unstructured  grids  either  with  a  subset  of  fine  grid  points  comprising  the  coarse  grids 
or  with  completely  independent  coarse  and  fine  meshes.  This  has  been  shown  to  be  successful  in  both 
for  inviscid  and  viscous  flow  computations  [105,  101,  129,  24].  Both  the  approaches  outlined  above 
share  a  common  problem,  that  of  generating  coarse  grid  levels.  For  complex  geometries,  especially  in 
three  dimensions,  generating  coarse  grids  that  faithfully  represent  the  complex  geometries  can  become 
a  difficult  proposition.  The  requirement  to  generate  not  just  one  grid  but  multiple  grids  that  preserve 
the  geometry  places  too  much  of  a  burden  on  a  user.  The  third  approach  circumvents  this  problem 
by  generating  coarse  grids  through  agglomeration  or  fusing  of  fine  grid  control  volumes,  resulting  in 
polyhedral  coarse  grid  control  volumes.  This  method  was  developed  in  [87]  for  cell-vertex  schemes  and 
in  [158]  for  cell-centered  schemes.  This  method  has  been  further  refined  to  deal  with  inviscid  and  viscous 
flows  past  complex  configurations  in  both  two  and  three  dimensions  [85,  177,  108,  109].  In  Figures  5  (a-e) 
the  results  from  using  the  agglomeration  multigrid  strategy  are  presented.  Figure  5a  depicts  the  surface 
grid  employed  to  compute  inviscid  flow  about  a  low-wing  transport  configuration  [177].  The  mesh 
contains  804,056  vertices  and  approximately  4.5  million  tetrahedra.  Figure  5b  displays  the  computed 
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surface  Mach  contours  for  transonic  inviscid  flow  using  a  seven-level  agglomeration  multigrid  strategy. 
The  freestream  conditions  for  this  case  are  M^o  =  0.77  and  1.116°  incidence.  Figure  5c  shows  the 
surface  mesh  employed  for  computing  turbulent  viscous  flow  over  a  partial  span-flap  wing  experimental 
configuration  [109].  The  fine  grid  contains  nearly  2.-3  million  vertices  and  nearly  13.6  million  tetrahedra. 
The  freestream  conditions  for  this  case  are  Moo  =  0.2,  incidence  of  10°,  and  Reynolds  number  of  2 
million,  corresponding  to  an  approach  condition.  The  solution  obtained  with  a  six-level  agglomeration 
multigrid  strategy  is  shown  qualitatively  in  Figure  5d  in  the  form  of  Mach  contours  on  the  wind-tunnel 
wall  and  density  contours  on  the  wing  surface.  Figure  5e  shows  the  convergence  histories  for  these  two 
cases.  The  deterioration  in  convergence  of  the  Navier-Stokes  case  is  mainly  due  to  the  use  of  stretched 
meshes  and  indicates  that  the  efficiency  of  Navier-Stokes  solvers  is  far  from  satisfactory. 


6.1.2  Implicit  schemes 

Implicit  schemes  for  the  compressible  Euler  and  Navier-Stokes  equations  have  been  developed  in  order 
to  accelerate  the  convergence  to  steady  state.  If  the  time  derivative  in  Eqn.  (7)  is  replaced  by: 

dW  -  IT” 

^ 

an  explicit  scheme  is  obtained  by  evaluating  R{W)  at  time  level  n.  An  implicit  scheme  is  obtained  by- 
evaluating  R{W)  at  level  n  +  1.  In  the  latter  case,  linearizing  R  about  time  level  n,  we  obtain 


ait;- =  (iF"+^  -  (10) 

Eqn.  (9)  represents  a  large  nonsymmetric  linear  system  of  equations  for  the  updates  of  the  vector 
of  unknowns  and  needs  to  be  solved  at  each  time  step.  As  At  tends  to  infinity,  the  method  reduces 
to  Newton’s  method.  Direct  solvers  have  been  used  solve  Eqn.  (9)  yielding  quadratic  convergence, 
but  these  entail  prohibitive  costs  and  memory  requirements  and  are  impractical  for  three  dimensional 
applications  [176,  157].  Typically,  due  to  storage  considerations,  only  a  lower  order  representation  of  the 
left-hand  side  is  employed  and  the  system  is  solved  by  iterative  means.  As  a  result  of  this  approximation 
Eqn.  (9)  can  never  approach  Newton’s  method  (with  its  associated  quadratic  convergence  property) 
due  to  the  mismatch  of  the  right  and  left  hand  side  operators. 

Thareja  et  al.  [163]  and  Hassan  et  al.  [58]  have  utilized  point-implicit  iterative  procedures.  Fezoui 
[45],  Batina  [17],  Anderson  et  al.  [5,  7]  and  Slack  et  al.  [157]  have  used  a  Gauss-Seidel  relaxation 
technique.  It  is  also  possible  to  use  more  sophisticated  techniques  for  the  solution  of  the  linear  system, 
such  as  GMRES  [149]  and  QMR  [46].  Preconditioning  of  the  matrices  is  critical  to  achieving  good 
convergence  with  these  methods.  Some  of  the  typical  preconditioners  are  block-diagonal,  symmetric 
successive  over- relaxation  (SSOR)  and  Incomplete  LU  factorization  (ILU)  preconditioners  [111].  GM¬ 
RES  with  diagonal  preconditioning  has  been  used  by  Shakib  et  al.  [151]  to  solve  the  linear  systems 
arising  out  of  a  finite  element  discretization  of  the  Euler  equations.  Slack  et  al.  [157]  and  Whitaker  [184] 
have  also  used  GMRES  with  diagonal  preconditioning  in  two  and  three-dimensional  applications.  Slack 
et  al.  [157]  have  observed  when  solving  the  two-dimensional  Euler  equations  that  the  preconditioned  it¬ 
erative  methods  perform  better  than  the  other  methods  as  the  number  of  elements  in  the  mesh  increases. 
Venkatakrishnan  and  Mavriplis  [178]  tested  a  family  of  implicit  schemes  for  solving  the  two-dimensional 
compressible  Navier-Stokes  equations  on  unstructured  meshes.  They  concluded  that  GMRES  with  ILU 
preconditioning  (GMRES/ILU)  was  superior  to  other  implicit  schemes  over  a  range  of  problems  sizes 
and  flow  conditions  and  was  competitive  in  terms  of  CPU  time  with  multigrid  methods.  Luo  et  al. 
[99]  used  GMRES/ILU  to  compute  two-  and  three-dimensional  flows.  A  drawback  of  all  the  methods 
based  on  linear  algebra  is  that  the  memory  requirements  become  severe  which  seriously  limits  the  sizes 
of  the  problems  that  can  be  solved,  especially  in  three  dimensions.  It  is  possible  to  implement  GMRES 


9 


in  a  matrix-free  form  [29]  where  the  matrix- vector  product  is  replaced  by  a  finite  difference  expression 
involving  residual  evaluations.  Another  advantage  of  the  matrix-free  approach  is  that  the  higher  order 
discretization  can  be  employed  for  the  left-hand  side  and  the  matrix- vector  product  [.30].  The  drawback 
of  this  approach  is  that  powerful  preconditioners  like  SSOR  and  ILU  cannot  be  used  since  they  require 
the  matrix  to  be  explicitly  available.  A  compromise  is  to  settle  for  a  matrix-free  form  of  GMRES  with 
diagonal  preconditioning  as  done  in  [75]. 

Another  implicit  method  that  has  been  investigated  is  based  on  the  use  of  linelets  [59].  The  key  idea 
is  to  cover  the  domain  with  a  set  of  lines  and  the  scheme  is  made  implicit  along  these  lines.  Thus  tri- 
diagonal  systems  need  to  be  solved  which  can  be  accomplished  efficiently.  A  drawback  of  this  method 
is  that  the  convergence  is  sensitive  to  the  orientation  of  the  linelets  with  respect  to  the  direction  of 

strongest  coupling. 

6.2  Unsteady  Problems 

While  the  solution  techniques  for  computing  steady  flows  have  evolved  to  a  high  degree,  those  for  dealing 
with  unsteady  flows  have  lagged  behind.  This  is  partly  due  to  fact  that  prediction  of  aerodynamic 
properties  in  steady  transonic  flow  (cruise  conditions)  has  been  the  driving  force,  with  the  unsteady 
effects  being  important  only  in  extreme  conditions,  such  as  flutter  and  stall.  It  is  anticipated,  however, 
that  unsteady  time-accurate  simulations  will  provide  the  next  great  challenge  for  CFD. 

Recall  that  after  discretization  in  space,  a  system  of  coupled  ODE’s  results  (Eqn  (6)).  If  we  assume 
that  the  system  is  decoupled  i.e.  M  =  I,  the  simplest  way  to  solve  the  system^  of  ODE’s  is  to  use  an 
explicit  scheme.  Low-storage  Runge-Kutta  methods  [67,  170]  and  various  predictor-corrector  schemes 
faU  under  the  category  of  expficit  schemes.  These  schemes  typically  require  only  simple  updates.  They 
place  restrictions  on  time  step  due  to  the  CFL  condition,  which  become  severe  in  the  case  of  semi¬ 
discrete  schemes  satisfying  monotonicity  principles.  However,  explicit  schemes  may  be  the  schemes  of 
choice  for  certain  unsteady  applications  when  the  time  scales  of  interest  are  small  or  more  precisely, 
that  they  are  comparable  to  the  spatial  scales.  The  grid  should  be  clustered  only  in  regions  of  interest; 
otherwise,  the  size  of  the  explicit  time  step  could  be  unnecessarily  small.  The  situation  can  be  improved 
by  the  use  of  time-step  sequencing  [83,  138],  where  different  ceUs  take  varying  number  of  local  time 
steps  to  get  to  a  particular  time  level. 

When  the  mass  matrix  is  present,  even  the  explicit  schemes  require  a  matrix  inversion.  Often,  either 
the  mass-matrix  is  lumped  for  convenience  or  at  best  a  few  Jacobi  iterations  are  carried  out  [124,  36]. 
Neither  approach  is  entirely  satisfactory,  especially  if  large  time  steps,  permitted  by  multi-stage  explicit 
schemes  and  implicit  schemes,  are  used.  Donea  [41]  realized  that  the  presence  of  the  mass  matrix  in 
his  formulations  wiU  make  it  unattractive  and  proposed  a  simple  two-pass  procedure.  While  he  proved 
that  the  formal  order  of  accuracy  remained  the  same,  he  also  showed  that  the  dissipation  and  dispersion 
properties  were  somewhat  compromised.  Recently,  it  has  been  observed  [179]  that  the  mass  matrix  can 
be  lumped  without  adverse  consequences  if  the  cell-vertex  scheme  possesses  only  second  or  lower  order 
accuracy. 

When  an  implicit  scheme  is  used  to  solve  for  unsteady  flows,  the  linear  system  Eqn. (9)  is  modified 
by  replacing  the  diagonal  matrix  V  by  the  matrix  MV  where  M  is  the  mass  matrix  and  I  is  the  volume 
of  the  cell.  However,  it  is  not  enough  to  solve  this  linear  system.  One  has  to  drive  the  unsteady  residual 
to  zero  or  at  least  to  truncation  error.  This  is  usuaUy  done  by  employing  inner  iterations  [139.,  136].  It 
is  the  role  of  these  inner  iterations  to  eliminate  errors  due  to  the  factorization  (if  any  is  carried  out), 
linearization,  and  errors  arising  from  employing  a  lower  order  approximation  on  the  implicit  side.  The 
number  of  inner  iterations  required  may  be  large  depending  on  the  flow  situation  and  the  size  of  the 

time  step  employed. 

Brennis  and  Eberle  [27]  and  Jameson  [68]  have  advocated  a  different  approach  for  deriving  an 
efficient  implicit  scheme  for  unsteady  flows.The  idea  is  to  define  an  unsteady  residual,  foUowdng  a 
backward  difference  approximation  to  the  time  derivative  and  then  use  either  a  relaxation  strategy  or  a 
multigrid  technique  to  drive  the  unsteady  residual  to  zero.  The  significant  advantage  of  this  approach 
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when  multigrid  is  used  to  solve  the  nonlinear  problem  is  that  it  incurs  no  storage  overheads  associated 
with  traditional  implicit  schemes,  and  is  particularly  attractive  for  unstructured  grid  computations  in 
three  dimensions.  It  allows  the  time  step  to  be  determined  solely  based  on  flow  physics.  This  method 
has  been  used  to  compute  two-  and  three-dimemsional  flows  over  airfoils  and  wings  [68,  112,  4]  using 
structured  grids.  Vassberg  [172]  has  applied  this  method  to  compute  flow  solutions  over  oscillating 
airfoils  using  unstructured  grids.  This  approach  has  also  been  used  in  [179]  in  conjunction  with  the 
agglomeration  multigrid  strategy  to  solve  two-dimensional  unsteady  flows  where  the  inversion  of  the 
mass  matrix  is  accompUshed  indirectly  during  the  multigrid  procedure. 


7  Adaptive  grids 

Solution  adapted  grids  are  increasingly  being  used  to  compute  complex  unsteady  and  steady  flows. 
There  are  three  distinct  ways  the  grid  can  be  adapted  to  the  solution.  These  are  r-refinement,  h- 
reflnement  and  p-reflnement.  In  r-refinement,  the  nodes  are  simply  redistributed  so  that  regions  of 
importance  are  better  resolved.  In  h-refinement  or  mesh-enrichment,  the  cells  are  locally  subdivided 
or  merged  or  in  some  instances,  a  complete  remeshing  is  done  to  reduce  the  grid  spacing  in  regions  of 
interest.  In  p-refinement,  the  degree  of  the  basis  function  is  adjusted  locally  to  match  the  variation 
in  solution.  For  a  survey  on  adaptive  mesh  refinement  techniques,  the  reader  is  referred  to  the  review 
article  by  Powell  et  al.  [1.35]. 

R-refinement  is  probably  the  simplest  in  concept,  but  is  burdened  with  practical  difficulties  in  multi¬ 
dimensions  especially  when  dealing  with  highly  stretched  grids.  For  a  survey  of  mesh  point  movement 
methods,  see  the  review  article  by  Hawken  et  al.  [60].  The  difficulties  include  skewness,  crossing  of 
lines,  arbitrarily  small  cell  volumes  etc.  Some  progress  has  been  made  in  dealing  with  these  issues 
at  least  for  inviscid  flows,  where  the  cell-aspect  ratios  are  not  too  severe  [122].  The  advantage  of  r- 
refinement  is  that  if  a  valid  grid  results  from  it,  all  that  is  required  is  the  interpolation  of  variables  from 
the  old  to  the  new  grid.  This  could  be  done  in  a  conservative  manner  if  desired.  A  way  to  avoid  the 
interpolation,  which  introduces  errors  that  could  accumulate,  is  to  introduce  the  grid  movement  terms 
in  the  governing  equations  Eqn.  (1).  These  terms  need  to  be  discretized  carefully  so  that  freestream 
is  preserved.  The  Geometric  Conservation  Law  (GCL)  [165,  186]  formalizes  this  procedure.  Recently, 
r-refinement  has  been  used  to  great  advantage  with  Roe’s  upwind  scheme  [143]  to  obtain  “fitted”  shock 
resolution  by  Paraschivoiu  et  al.  [123],  Parpia  and  Parikh  [126],  and  van  Rosendale  [171]  by  aligning 
edges  of  the  triangulation  with  discontinuities.  The  “fitting”  is  done  in  a  shock-capturing  framework 
by  utilizing  that  property  of  Roe’s  scheme  which  allows  isolated  discontinuities  aligned  with  the  mesh 
to  be  captured  exactly. 

H-refinement  is  by  far  the  most  popular  means  of  adaptation  in  compressible  flows.  This  is  especially 
true  for  inviscid  flows  dominated  by  interactions  of  shock  waves  where  p-refinement  techniques  are  of 
limited  value.  In  the  Adaptive  Mesh  Refinement  (AMR)  approach,  h-refinement  is  employed  with  body- 
fitted  structured  grids  to  adapt  to  the  flow  solution  [20,  138].  Instead  of  body-fitted  grids,  Cartesian 
meshes  can  be  used  to  adapt  to  complex  geometries  in  addition  to  adapting  to  flow  features  [21,  37,  137]. 
Since  some  of  the  approaches,  e.g.  [37,  152]  do  not  make  use  of  any  structure  and  employ  indirect 
addressing,  these  methods  may  be  classified  as  unstructured.  The  AMR  framework  has  been  used 
to  compute  many  complex  inviscid  unsteady  flows  that  exhibit  disparate  spatial  and  temporal  scales. 
Unstructured  grids,  composed  of  triangles  or  tetrahedra,  are  also  particularly  suited  for  h-refinement.  In 
the  case  of  both  AMR  and  unstructured  grids,  the  regions  of  interest  are  first  identified  either  through 
a  combination  of  heuristic  criteria  such  as  density  gradients  (undivided)  or  through  estimation  of  the 
truncation  error.  In  the  case  of  AMR,  typically  cells  are  subdivided  in  some  fixed  ratio,  and  these  are 
sometimes  grouped  into  logically  rectangular  patches.  In  the  case  of  unstructured  grids,  the  regridding 
procedure  is  more  involved.  If  a  Delaunay  triangulation  is  used  as  the  initial  grid,  one  can  use  Watson’s 
incremental  algorithm  [183]  to  obtain  a  new  Delaunay  triangulation.  Other  procedures  such  as  edge¬ 
swapping  may  also  be  used  to  restore  the  Delaunay  tri angulation.  In  the  case  of  other  triangulations, 
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cells  are  subdivided  and  smoothed  using  a  few  passes  of  a  Laplacian-type  smoother.  Peraire  et  al.  [130] 
advocate  the  regeneration  of  the  grid  by  the  advancing  front  method  for  steady  state  applications  after 
identifying  regions  of  interest  by  error  estimation.  In  the  case  of  nodal  schemes  in  three  dimensions, 
when  edges  are  tagged  and  the  cells  are  to  be  refined  or  derefined  based  on  these  tagged  edges,  the 
possibilities  are  numerous  [79].  It  helps  to  use  some  rules  to  limit  the  cases  [35].  Several  investigators 
have  employed  multigrid  procedure  for  steady  state  computations  using  the  grids  generated  by  solution 
adaptation  [102,  127,  35,  26].  Figure  6a  shows  the  initial  surface  grid  for  an  isolated  engine  inlet. 
The  initial  grid  contained  3536  nodes.  After  three  adaptive  refinement  cycles,  the  final  mesh  shown  in 
Figure  6b  was  obtained.  The  final  mesh  contained  61,002  nodes.  This  computation  done  by  Connell 
and  Holmes  [35]  utilized  the  adaptive  multigrid  procedure.  When  adapting  grids  to  flow  features,  a 
problem  that  may  occur  is  that  certain  regions  are  overly  resolved  at  the  expense  of  other  regions, 
resulting  in  seemingly  sharp  but  incorrect  solutions.  Warren  et  al.  [182]  have  proposed  and  tested  a 
simple  modification  to  the  error  indicators  that  alleviates  this  problem. 

In  the  case  of  transient  problems,  adaptation  is  performed  much  more  frequently  and  therefore  the 
regridding  process  is  required  to  be  efficient.  While  this  is  easily  accomplished  efficiently  by  subdivision 
of  cells  in  the  AMR  framework,  the  problem  is  more  involved  in  the  case  of  unstructured  grids  [94,  141]. 
The  problem  of  flow  past  bodies  in  relative  motion  has  also  been  addressed  in  the  literature  [93,  167, 
57.  82,  31].  Typically,  mesh  point  movement  and  efficient  mesh  restructuring  are  employed  to  obtain 
valid,  good-quality  grids  about  the  moving  bodies.  Figures  7(a-c)  show  the  results  from  the  simulation 
of  store  separation  from  F-117.  The  freestream  conditions  are  Moo  =  0-84  and  incidence  of  2.73°. 
The  projectiles  are  forcibly  ejected  at  different  velocities.  This  multi-body  problem  was  computed  by 
Baum  et  al.  using  the  Arbitrary  Euler  Lagrangian  adaptive  unstructured  methodology  developed  in 
[18].  Because  of  symmetry,  the  motion  of  only  two  of  the  projectiles  is  simulated.  Figures  7a  and  7b 
show  the  Mach  contours  from  two  different  perspectives  and  Figure  7c  displays  the  pressure  contours 
at  a  particular  instant  of  motion.  The  complex  shock  patterns  emanating  from  the  projectiles  and  their 
signatures  on  the  plane  of  symmetry  may  be  observed  in  Figure  7c. 


8  Higher  order  methods 

The  use  of  methods  possessing  higher  than  second  order  accuracy  for  solving  the  compressible  Navier- 
Stokes  equations  on  unstructured  grids  is  not  yet  commonplace  and  is  a  topic  of  active  research.  The 
finite  element  method  achieves  higher  order  accuracy  by  the  use  of  higher  order  elements.  Bey  and 
Oden  [22]  have  used  a  discontinuous  Galerkin  method  to  obtain  up  to  4th  order  accuracy  for  smooth 
fiow^s  using  structured  adaptive  grids.  In  the  finite  volume  community,  Barth  and  Fredrickson  [14] 
derived  general  conditions  for  a  scheme  to  be  higher  order  accurate  that  involve  the  reconstruction 
of  variables  satisfying  the  properties  of  conservation  of  mean,  ^-exactness  and  compact  support,  k- 
exactness  refers  to  the  property  of  being  able  to  reconstruct  exactly  polynomials  of  degree  <  k.  The 
key  idea  is  to  extend  the  support  of  the  scheme  to  enable  the  coefficients  in  the  polynomial  to  be 
determined.  They  also  proposed  a  minimum-energy  reconstruction  that  used  a  larger  support  and 
utilized  a  least-squares  method  to  evaluate  the  coefficients.  By  combining  this  reconstruction  procedure 
with  higher  order  quadrature,  they  demonstrated  the  effectiveness  of  both  h-  and  p-refinements  in 
reducing  the  errors  in  smooth  flows.  More  recently,  Chakravarthy  et  al.  [31]  have  proposed  a  similar 
approach  to  achieving  higher  order  accuracy.  Equivalent  to  requiring  /^-exactness,  they  propose  that 
the  mean  be  conserved  in  a  neighborhood  of  the  point  in  question.  Mitchell  [115]  has  proposed  a  quasi¬ 
quadratic  reconstruction  procedure  that  appears  to  yield  better  results  than  those  obtained  with  linear 
reconstruction.  Harten  and  Chakravarthy  [55]  have  proposed  a  framework  for  applying  Essentially 
Nonosdllatory  (ENO)  schemes  [56]  on  unstructured  meshes.  They  present  two  techniques  for  adapting 
the  stencils  during  the  reconstruction  procedure,  which  also  give  preference  to  a  centered  stencil  in 
smooth  regions.  The  first  technique  considers  several  candidate  stencils  and  picks  the  one  that  represents 
the  function  the  smoothest.  The  second  technique  is  a  hierarhical  one  that  incrementally  augments  the 
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stencil  as  derivatives  are  required  during  the  reconstruction  procedure.  Abgrall  [2]  has  also  proposed  a 
hierarchical  ENO  scheme  and  has  tested  it  by  computing  a  variety  of  inviscid  compressible  flows. 

A  problem  with  the  reconstruction  methods  described  above  is  that  since  the  support  grows  with  the 
order  of  the  scheme,  computational  cost  and  boundary  treatment  become  significant.  Halt  and  Agarwal 
[.5.3]  have  proposed  and  tested  two  methods  for  realizing  higher  order  accuracy  while  keeping  the  support 
compact.  These  techniques  are  based  on  deriving  equations  for  either  the  derivatives  or  the  moments 
of  the  governing  equations.  They  demonstrate  up  to  4th  order  accuracy  with  both  these  methods  for 
smooth  flows.  They  conclude  that  the  method  based  on  moments  is  more  suitable  for  dealing  with  shocks 
and  boundary  conditions.  Barth  [13]  has  also  proposed  a  different  quadratic  reconstruction  procedure. 
Similar  to  a  finite  element  procedure,  additional  degrees  of  freedom  are  placed  at  the  mid-points  of  edges 
in  a  triangulation.  The  gradient  and  Hessian  required  for  quadratic  reconstruction  are  again  obtained 
by  a  least-squares  procedure.  Figures  8(a-d),  provided  by  Barth  [13],  illustrate  the  results  obtained  with 
various  discretizations  for  the  circular  advection  problem: 


Ut  -f-  {^yu^x 

Discontinuous  initial  data  is  specified  along  an  interior  cut  line.  The  exact  solution  is  a  solid  body 
rotation  of  the  cut-line  data.  The  triangular  grid  is  displayed  in  Figure  8a.  Solutions  obtained  with 
the  finite  volume  scheme  employing  piecewise-constant,  piecewise-hnear  and  piecewise-quadratic  recon¬ 
structions  are  depicted  in  Figures  8  (b-d).  The  significant  improvement  with  higher  order  reconstruction 
is  apparent.  Barth  [13]  also  states  that  quadratic  reconstruction  is  about  7  times  as  expensive  as  linear 
reconstruction;  of  this  a  factor  of  4  arises  from  the  introduction  of  grid  points  at  the  mid-points  of  the 
edges.  This  “h-refinement”  contributes  in  part  to  the  improvement  realized  with  quadratic  reconstruc¬ 
tion. 

9  Hybrid  discretizations 

While  triangulations  in  two  and  three  dimensions  offer  flexibility  in  adapting  to  complex  geometries, 
it  is  not  clear  if  they  are  needed  near  the  wall  regions.  Steger  [161]  when  introducing  the  “thin- 
layer”  approximation  in  1978  stated,  “One  does  not  generally  have  sufficient  computer  to  resolve  the 
viscous  terms  except  in  a  thin-layer  near  the  body”.  This  statement  holds  true  even  today  for  practical 
aerodynamic  computations.  Thus  highly  stretched  triangulations  are  required  in  the  near-wall  regions. 
They  can  be  generated  by  employing  a  mapping  procedure  [10.3]  where  a  Delaunay  triangulation  is 
performed  in  the  mapped  plane  which  is  then  mapped  back  to  the  physical  space  to  yield  the  desired 
stretched  triangulation.  This  method  has  also  been  used  for  adaptation  in  viscous  flows  [181].  Aftosmis 
et  al.  [3]  have  addressed  the  important  issue  of  element  shapes  by  examining  the  accuracy  obtained 
on  test  problems  by  using  various  triangular  element  shapes  and  quadrilaterals.  They  conclude  that 
using  stretched  triangulations  near  the  body  does  not  yield  any  improvement  in  solution  accuracy 
over  using  quadrilaterals,  and  is  therefore  inefficient  because  of  the  extra  edges  in  the  triangulation. 
They  recommend  a  method  that  removes  unnecessary  edges  (diagonals)  from  boundary  layer  regions. 
Rather  than  adopt  this  two-step  procedure,  where  the  triangular  grid  is  first  generated  followed  by  the 
ehmination  of  unnecessary  edges,  the  hybrid  method  ehminates  the  need  for  stretched  triangulations 
all  together.  This  method  makes  use  of  a  structured  or  a  semi-structured  grid  in  the  near-wall  regions. 
In  two  dimensions,  a  structured  body-fitted  grid  is  used  near  the  body  whereas  in  three  dimensions 
a  prismatic  grid  is  used  since  the  surface  grid  is  triangular.  Nakahashi  et  al.  [117,  118],  KaUinderis 
et  al.  [80,  128],  Connell  and  Braaten  [-34]  have  advocated  this  approach.  The  generation  of  prismatic 
elements  requires  careful  marching  away  from  the  surface  so  that  lines  do  not  cross.  Melton  et  al.  [11-3] 
have  used  a  prismatic  grid  with  a  background  Cartesian  mesh  to  compute  three-dimensional  flows.  The 
structure  afforded  could  be  effectively  used  to  tailor  line-implicit  algorithms,  although  this  requires  that 
the  various  elements  be  processed  separately.  Ideally,  one  data  structure  (edge-based)  should  be  able  to 
process  aU  elements  regardless  of  their  shapes  for  an  automatic  procedure.  Figure  9a  shows  the  surface 
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grid  for  an  adapted  hybrid  prismatic-tetrahedral  grid  about  a  hemisphere;  The  initial  grid  contained 
3000  triangular  faces  on  the  wall  and  20  layers  of  prisms  and  70,000  tetrahedra.  The  freestram  conditions 
for  supersonic  laminar  flow  are  Moo  =  1-4  and  Re  =  1000.  After  the  dual  adaptation  procedure,  where 
the  prisms  are  directionally  refined  with  the  normal  resolution  unaltered  and  the  tetrahedra  are  divided 
into  either  2,4,  or  8  tetrahedra,  the  final  grid  contains  7000  triangular  faces  and  20  layers  of  prisms  and 
225,000  tetrahedra.  The  signature  of  the  directional  refinement  of  the  prisms  may  be  observed  on  the 
hemispherical  surface.  Figure  9b  displays  the  surface  Mach  contours  on  the  symmetry  plane  for  this 
flow.  These  figures  have  been  provided  by  Parthasarathy  et  al.  [128].  The  explicit  solver  also  makes 
use  of  an  adaptive  multigrid  procedure. 

10  Parallel  computing 

Computational  fluid  dynamics  as  its  name  implies  is  inevitably  linked  to  computing  issues.  Among 
these  are  processing  power,  memory  technology,  networking  and  accessibility.  Ability  to  compute  the 
solutions  to  problems  in  finite  time  always  being  the  goal,  CFD  has  benefited  immensely  from  the  rev¬ 
olution  that  has  taken  place  in  the  last  15  years  in  these  areas.  Vector  supercomputers  have  provided 
much  of  the  computing  power  that  has  been  harnessed  to  compute  complex  three-dimensional  flows. 
It  is  anticipated  that  distributed-memory  parallel  computers  will  offer  the  next  cost-effective  leap  for¬ 
ward  in  terms  of  computing  power.  Much  research  has  been  carried  out  in  the  area  of  unstructured 
grid  computations  on  parallel  computers.  Some  of  the  issues  that  arise  are  partitioning  of  the  grid, 
message  patterns,  data  structures  and  design  of  parallel  algorithms.  Partitioning  of  unstructured  grids 
for  parallel  computing  has  been  investigated  by  a  number  of  researchers.  The  methods  can  be  broadly 
classified  into  geometry-based  [23,  114,  44]  and  graph-based  [133,  156]  algorithms.  It  appears  that  the 
graph-based  algorithms,  in  particular,  the  spectral  bisection  technique  of  Pothen  et  al.  [133],  while 
being  computationally  more  expensive,  yield  much  better  partitions.  Unstructured  grid  flow  solvers 
have  been  implemented  on  various  machines  [180,  107,  75,  43,  140].  These  studies  have  shown  that 
good  performance  may  be  obtained  by  paying  careful  attention  to  the  issues  listed  above.  Regarding 
the  flow  solver,  explicit  schemes  and  point  Jacobi  implicit  schemes  possess  almost  complete  parallelism, 
except  for  communication  at  the  inter-processor  boundaries.  Multigrid  methods  appear  to  have  ade¬ 
quate  parallelism  for  coarse-grained  parallel  computers;  impressive  performances  have  been  reported  by 
Mavriplis  et  al.  [107].  Johan  et  al  [75],  Ramamurthi  et  al.  [140]  and  Venkatakrishnan  [175]  have  shown 
that  implicit  schemes  can  be  carefully  designed  to  yield  good  performance  when  solving  unstructured 
grid  problems  on  parallel  computers.  While  a  matrix-free  GMRES  with  nodal  block  preconditioner 
that  has  almost  complete  parallelism  except  for  the  communication  at  the  inter-processor  boundaries 
is  used  in  [75],  methods  that  are  implicit  within  processors  but  explicit  across  processors  are  used  in 
[140,  175].  Thus,  as  the  number  of  processors  increases,  the  latter  methods  exhibit  a  degradation  in 
convergence.  In  [175]  it  is  demonstrated  that  convergence  rate  of  such  a  method  can  be  improved  by 
using  in  addition  a  global  coarse  grid  approximation,  where  each  processor  is  represented  by  one  grid 
point.  In  general,  it  appears  that  the  best  parallel  algorithm  is  also  the  best  serial  one,  with  some 
possible  loss  in  convergence  in  the  case  of  imphcit  schemes  in  a  multi-processor  environment. 


11  Conclusions 

In  this  paper,  unstructured  grid  flow  solvers  for  the  compressible  Navier-Stokes  equations  have  been 
surveyed.  Significant  progress  has  been  made  in  the  areas  of  spatial  and  temporal  discretization,  adaptive 
and  parallel  algorithms.  Newly  developed  field  equation  turbulence  models  seem  to  mesh  nicely  with 
the  unstructured  grid  framework.  Based  on  these  developments,  it  is  safe  to  say  that  unstructured 
grid  technology  is  almost  on  par  with  structured  grid  technology,  although  encumbered  with  additional 
memory  and  computational  costs.  This  overhead  has  to  be  balanced  with  the  ability  to  compute  flows 
over  complex  geometries  and  the  ease  of  adaptation. 
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The  areas  that  require  further  research  are  thus  the  same  as  in  structured  grids.  These  include  better 
and  faster  implicit /multigrid  flow  solvers  that  are  insensitive  to  cell  aspect  ratios  and  grid  stretching, 
improved  higher  order  discretization  techniques,  better  a  priori  or  a  posteriori  error  estimates  and 
parallel  algorithms.  The  status  of  unsteady  flow  solvers  is  far  from  satisfactory  either  for  structured  or 
unstructured  grids.  The  use  of  hybrid  meshes  needs  to  be  explored  further  to  utilize  data  structures 
that  can  handle  any  type  of  tessellation.  As  unsteady  flows  and  aeroacoustics  are  emerging  as  areas 
of  interest,  nondissipative  schemes  that  also  minimize  dispersion  have  to  be  designed  to  be  applicable 
to  unstructured  grids,  which  are  nonuniform  in  general.  In  this  respect  as  well  as  in  the  areas  of  error 
estimation  and  higher  order  schemes,  there  should  be  more  interaction  between  finite  volume  and  finite 
element  communities. 
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Figure  1.  Mach  contours  for  inviscid  transonic  flow  over  Boeing  747-200  using  an  unstructured  tetra¬ 
hedral  grid  {Moo  =  0.8,  a  =  2.7.3°)  [72]. 


Figure  2.  Laminar  flow  over  a  delta  wing  (M,^  =  0.3,  a  =  20.5°,  Re  =  0.9  x  10®)  [48].  (a)  Surface  grid 
and  “oil-flow'’  pattern  (b)  Pressure  profiles  at  four  chord  stations. 


station  A 


Station  A 


Figure  4.  Compressible  turbulent  flow  over  a  high-lift  configuration  =  0.2,  a  =  16°)  [6].  Surface 

pressure  distributions  and  velocity  profiles  at  various  stations. 


Fig.  5e. 

Figure  5.  Computations  done  using  the  agglomeration  multigrid  strategy  [177,  109].  (a)  Surface  grid  for 
the  low  wing  transport  (LWT)  case,  (b)  Mach  contours  on  the  surface  for  inviscid  transonic  flow  over 
the  LWT.  (c)  Surface  grid  for  the  partial  span-flap  experimental  configuration,  (d)  Mach  contours  on 
the  wind-tunnel  wall  and  density  contours  on  the  surface  for  turbulent  flow,  (e)  Convergence  histories 
for  the  inviscid  and  the  viscous  cases. 
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(a) 


(b) 


Figure  Y.  Store  separation  from  F-llY  [18].  Solution  at  a  particular  instant  in  time,  (a)  Mach  contours 
from  outside,  (b)  Mach  contours  from  inside  the  bay.  (c)  Pressure  contours  illustrating  the  shock  wave 
patterns  about  the  projectiles. 
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(c) 


(d) 


Figure  8.  Comparison  of  various  discretization  methods  for  the  circular  advection  problem  [13].  (a) 
Grid  used  for  the  two-dimensional  advection  equation,  (b)  Piecewise-constant  finite- volume  solution, 
(c)  Piecewise-linear  finite- volume  solution,  (d)  Piecewise-quadratic  finite- volume  solution. 
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